Excess Vibrational Modes and the Boson Peak in Model Glasses 
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The excess low-frequency normal modes for two widely-used models of glasses were studied at 
zero temperature. The onset frequencies for the anomalous modes for both systems agree well 
with predictions of a variational argument, which is based on analyzing the vibrational energy 
originating from the excess contacts per particle over the minimum number needed for mechanical 
stability. Even though both glasses studied have a high coordination number, most of the additional 
contacts can be considered to be weak. 
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Our understanding of liquids is based on the idea that 
hquid structure is largely determined by strong short- 
ranged repulsions and that the longer-ranged attractions 
can be treated as a perturbation Similar considera- 
tions were used to study jamming at zero temperature as 
a function of density [1, 0] . When only finite-ranged 
repulsions are included, so that non-overlapping parti- 
cles do not interact, there is a sharp jamming transi- 
tion at (j)c. In such systems, the low-frequency normal 
modes of vibration in the marginally jammed state are 
fundamentally different from the long-wavelength plane 
waves expected from elasticity theory 0, The un- 
usual nature of the modes is reflected in the density of 
vibrational states, D{uj), versus frequency, w 0, [!]■ At 
the transition, D{u!) has a plateau extending to w = 0. 
This is very different from the expected Debye scaling 
D{uj) cx ui'^^^ in d dimensions [J, which is normally ob- 
served in solids. At densities above 0c, the plateau no 
longer extends to w = but terminates at a frequency 
uj*. The dramatic rise in D{uj) at w* corresponds to the 
onset of anomalous modes. Similar excess modes, whose 
origins are still debated [1, i, 0, i, H [H [HI, [ll| , are ob- 
served in glasses. The anomalous modes in glasses giv e 
rise to what is known as the boson peak [l^, [3, [3, [3] , 
and are believed to be responsible for many characteristic 
low-temperature phenomena (Tt} . 

In this paper, we demonstrate that the theoretical 
framework [H, [l^, used to explain the modes in 
marginally jammed solids at zero temperature can be 
extended to the anomalous modes in two more realis- 
tic models of "repulsive glasses," by which we mean sys- 
tems that undergo glass transitions due to the repulsive 
part of their potentials as the temperature is lowered. 
The zero-temperature jamming transition of spheres with 
finite-ranged repulsions, which we will call Point J to dis- 
tinguish from other jamming transitions, coincides with 
the density at which the system has just enough contacts 
to satisfy the constraints of mechanical equilibrium 0. 
This is called the isostatic condition. The average coor- 
dination number z {i.e. the average number of particles 
with which a given particle interacts) needed for mechan- 
ical stability is Zc = 2d, where d is the spatial dimension 



[2ll |22| . By compressing the system or increasing the 
range of interaction, we increase z. These extra contacts 
suppress anomalous modes at low frequencies [18. J^, 20i] . 
We will show that in systems with long-ranged inter- 
actions, there is a well-defined division between a rel- 
atively few strong repulsive interactions (stiff contacts) 
and the more numerous weaker interactions (including 
attractions), which can be treated as a correction. 

We will study the onset frequency of the anomalous 
modes, uj\ in two widely-used models of glasses and the 
glass transition and cornpare the simulation results with 
theoretical predictions [2^. We use a mixture of 800 A 
and 200 B spheres with equal mass m interac ting in three 
dimensions via the Lennard- Jones potential [23[ : 
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where rij is the separation between particles i and j 
and eij and ct^ depend on the type of particles un- 
der consideration: sab — l-Se^A, ess = O.Se^A and 
o'AB = 0.8aAA, and asB = O.SSaAA. The poten- 
tial is cut off at Tij ~ 2.h(Tij and shifted to satisfy 
V{2.^(jij) — V'{2.^aij) = 0. For these systems, the den- 
sity is given in units of p = Nu\j^/ and e^A, ctaa and 
m are set to unity. We restrict the densities to lie above 
p « 1.2, where the pressure is positive. When this is not 
the case, there can be low-frequency modes arising from 
rather different physics [23 |. 

For comparison we also study a system with purely 
repulsive interactions where (j)c exists. We simulate a 
mixture of 500 A and 500 B spheres with as — l.icTA 
and equal mass m, interacting via the repulsive Lennard- 
Jones potential [l| 
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where e is the characteristic energy, and ct^ = {ai+(7j)/2. 
For these systems, we characterize density with the pack- 
ing fraction (j> = TT/6{NAcr^ + NBcr%)/L^. We set e = 1, 



2 



m = 1, and cta = 1- The simulations were all carried 
out in a cubic box with periodic boundary conditions. 
We study zero-temperature (T — 0) configurations which 
were obtained by quenching initially random (T = oo) 
configurations of particles to their local energy minima 
at T = using conjugate gradient energy minimization 

In the harmonic expansion, the energy of particle dis- 
placements SRi from their equilibrium positions is: 



SE = 



V'{j 



(3) 



where V and V" are the first and second derivatives of 
the pair potential Vir) with respect to r, and rij — "^ij^ij 
is the equilibrium separation vector between particles i 
and j. Here, 5Rij = SRi~5Rj, and SR:^ is the projection 

of SRij on the plane perpendicular to f y . We diagonalize 
the dynamical matrix [J] to obtain normal modes |n), and 
their corresponding eigenvalues The normal-mode 
frequencies are = where the index n runs from 1 
to Nd. From these frequencies, we calculate the density of 
vibrational states per particle, D{uj). It is instructive to 
also calculate D^^\uj), obtained by neglecting the second 
term in Eq. e. the stress term). This corresponds to 
replacing V{r) with unstretched springs, each chosen to 
have the same stiffness V"{r) as for the original potential. 

In Fig.[lja) and (b), we show two curves, D{uj) (heavy 
curves) and D'-^^uj) (fight curves), corresponding to the 
stressed and unstressed cases, at each density. Fig. [T](a) 
shows results for the repulsive Lennard-Jones systems. 
Close to the unjamming transition, these systems are 
approximately equivalent to the repulsive harmonic sys- 
tems studied previously Just above the transition, 
D{uj) has a plateau down to lu = 0. Upon compression 
the onset of these anomalous modes shifts up to w*, be- 
low which D{uj) decreases to zero. Fig. [Ub) shows the 
results for the system with Lennard-Jones interactions. 
Because there are attractive interactions, the jamming 
transition lies inside the liquid- vapor spinodal [2, and 
is inaccessible. Thus, the plateau in the density of states 
never extends to uj — 0. 

It is important to note that D^°^(w) falls much more 
sharply than D{lo) as uj decreases below uj* , but that 
D{uj) and D^^^{u>) are nearly indistinguishable above 
LU* . Here, uj* corresponds to the onset of the anoma- 
lous modes in unstressed systems. At high densities, 
the sharp peaks in D^^^{uj) below uj* arise from linear 
combinations of nearly-degenerate long- wavelength plane 
waves. (The lowest peak has plane waves with wavevec- 
tor, ki = 2tt/L where L is the box size, while the second 
peak has k2 — 2-\/27r/i as one would expect for the two 
lowest frequency modes in an elastic solid.) Thus, the 
onset of anomalous modes lies above these peaks. In the 
infinite system-size limit, the peaks should smooth out to 



yield the normal scaling: D{u>) 



UJ 



in d dimensions. 
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FIG. 1: Density of vibrational states per particle of 3D glasses 
of N = 1000 particles. We plot D{u) for stressed sys- 
tems (heavy curves) and D'°-'(tj) for unstressed systems (light 
curves) interacting via (a) the repulsive Lennard-Jones poten- 
tial at A(l> = lO"*^ (solid), 10"^ (dotted), 0.1 (dot-dashed), 
and 0.2 (dashed); and (b) the Lennard-Jones potential at 
p = 1.2 (sohd), 1.4 (dotted) and 1.6 (dashed). 



eral, extra contacts increase the frequency of the lowest- 
frequency anomalous modes in two ways: (i) they can 
increase the energy cost of a mode by adding extra nodes 
so that some bonds are unstretched during an oscillation, 
and (ii) they can leave the number of nodes fixed but in- 
stead increase the average restoring force (and therefore 
the energy) for the normal-mode displacement. 

A variational argument calculates an upper bound for 
the energy of a normal mode by minimizing the energy 
with respect to these two contributions. The first term 
in Eq. [3] indicates that a good trial function would have 
nodes, i. e., small values of {SRij ■fij)'^, where V" is large. 
Thus the ziN/2 contacts with the highest values of V" 
introduce nodes in the trial function. The remaining (z — 
zi)N/2 contacts increase the energy of the trial mode by 
increasing the restoring force for displacements according 
to the first term in Eq. [3] |l20[. We rewrite Eq. [3]as 



5E^°^ ^ SEi + 6E2, 
5E = 5E^°'^ +SE3, 



(4) 



We now recapitulate the theoretical ideas that address 
the properties of high-coordination systems [20j . In gen- 



where SEi, SE2 and SE^ represent the energy costs asso- 
ciated respectively with the ziN/2 stiffest contacts, the 
remaining (z — zi)N /2 weak contacts, and the contribu- 
tion of the stress term (the second term in Eq. [31 see 
Ref. [ii]). Note that 5E''°^ is the energy cost of a mode 
for the unstressed system. 

We now construct approximate expressions for these 
contributions. Ref. argues that « Aiki{zi — ZcY, 
where fci = (y")i is the average over the ZiN/2 contacts 
with the highest values of V" . This is consistent with 
earlier simulation results for harmonic springs for zi — 
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•2c < 3 [19|. While ki varies strongly with density and 
potential, we expect Ai to depend only weakly on these 
quantities [26l |. To obtain Ai, we compare to simulations 
of unstressed systems with springs of equal stiffness. The 
precise value of Ai depends on which point in the density 
of states we choose to represent the onset of anomalous 
modes. In the following analysis, we use the value uj^ 
where Z?(*''(ci;) reaches 0.25 of its maximum height, which 
fixes Ai — 0.018. We would have obtained a somewhat 
different constant if we had compared with data at a 
different point in the rise, but the difference would be 
small because D^^\uj) rises abruptly. 

To estimate 6E2, we assume that for the normalized 
trial mode, the displacements of i and j are uncorrelated 
with each other if i and j are connected by weak con- 
tacts. Then {{SRij ■ fj,-)^) ^ 2/Nd for an iV-particle sys- 
tem in d dimensions [20] ■ Thus SE2 ~ Y^^j ^" (''^v ) ' 
where the sum ^™ runs only over pairs of particles i 
and j connected by the {z — zi)N/2 weakest contacts. 
We have shown numerically that this approximation is 
reasonable for a system with particles connected by har- 
monic springs at their equilibrium lengths with two very 
different stiffnesses. For such a system, the stiff springs 
contribute only to 5Ei , the weak springs contribute only 
to (5i?2- Our expressions for 5E2 and 5Ei were thus ver- 
ified cleanly. 

We estimate 5E^ as follows. For uncorrelated displace- 
ments between particles i and j, {{5R^Y) « 2{d—l)/dN, 

leading to SE3 = A3S, where S = l/Nj^tj -^^tt^' ^^'^ 
A3 « (d— l)/d is of order unity. 

Note that 6E3 does not depend on zi, so it does not 
affect the energy minimization with respect to zi or zi 
itself. In order to compute zi , we thus compute the total 
energy cost in absence of stress: 
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5E^"'> = 0.018fci(zi - zcf 
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where ki = {V")i is the average over the ziN/2 contacts 
with the highest values of V" , and sums over the 
pairs of particles connected by the {z — zi)N/2 contacts 
with the smallest values of V" . By minimizing SE^°^ with 
respect to zi, we obtain an upper bound on the energy 
and therefore the frequency, u* — \/SE^^, of the lowest- 
frequency anomalous modes in the unstressed systems. 

In Fig. [21 we compare this theoretical prediction for u* 
with the onset of anomalous modes, uj\ from simulations 
in the unstressed system (determined by where D^^^ {uj) 
reaches 0.25 of its maximum value). The agreement is 
excellent everywhere, with no adjustable parameters. 

Our results show that even though a typical repulsive 
amorphous solid may have a high coordination number 
and therefore appear to be far from the unjamming tran- 
sition, most of the contacts are weak. The number of stijf 
contacts is only slightly in excess of the minimum needed 
for mechanical stability. This is shown in Fig. [21 (solid tri- 
angles) for both potentials at all densities studied. Even 




1.2 -(b) 
S 0.8 - 
® 0.4 - 




0.6 



to 



1.1 



1.3 



1.5 



1.7 



0.2 




FIG. 2: Left axis: The characteristic frequency oj^ (open 
symbols) calculated from numerical simulations and the cor- 
responding theoretical predictions uj* (crosses). Lines are to 
guide the eye and connect the theoretical points. Right axis: 
the fractional deviation from isostaticity SZi = {zi — Zc)/zc 
(closed symbols), as functions of (a) distance from the unjam- 
ming transition, — (j>c, for repulsive Lennard Jones mixtures 
and (b) density p for Lennard- Jones mixtures. 



for the Lennard-Jones system, where the power-law tail 
of the interaction leads to a divergent total coordination, 
we find that {zi — Zc)/zc < 0.6 at all densities. Thus 
{zi — Zc)/zc is a small parameter. 

While the results of Fig. [21 are for unstressed systems, 
real glasses have nonzero stresses. To estimate the onset 
of anomalous modes for systems with stress, we must 
calculate the total energy cost of a mode, 6E = (5£''"' -I- 
5E3. Fig. [21 shows that for our systems, SE'^^^ « —SE3, 
where JE''"' is evaluated at its minimum with respect to 
zi for both potentials at different densities. Thus, there is 
a near cancellation of two large terms leading to a small 
value of SE (with a large uncertainty) and therefore a 
very low onset frequency for the anomalous modes in the 
stressed systems. This is consistent with the finding that 



the expected scaling of £'(w) 



for plane- wave normal 



modes is eclipsed by an approximately linear frequency 
dependence at small uj for both potentials. The near- 
cancellation may be a result of the history of how the 
system was prepared [Toj . 

Our results provide a plausible explanation for the ori- 
gin of the excess vibrational modes of the boson peak 
in repulsive glasses. For two commonly-studied mod- 
els, we have shown that the boson peak derives from 
the same low-frequency anomalous modes that arise at 
Point J for systems with finite-ranged repulsions. Sev- 
eral theories have been advanced previously for the bo- 



son peak I, i, [3, i, i, M 
as well as disordered ones 




Fractal systems Q 
can exhibit excess 
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FIG. 3: The energy cost of the lowest frequency anoma- 
lous mode for an unstressed system, as a function of 
S = i/N '^^.V' {rij)/rij for repulsive Lennard- Jones mix- 
tures (solid circles) and Lennard-Jones mixtures (open trian- 
gles). The contribution of the stress term to the energy cost 
of a mode is SE^ — A^S, where A3 is of order unity. The 
straight line fit corresponds to SE^'^^ = —1.65. 

vibrational modes, but glasses are typically not frac- 
tal and low-coordination crystals can also display excess 
modes [11] . Approaches by Phillips \27\ , Thorpe .2^ and 
Alexander [24I are also based on the idea that a mini- 
mum average coordination number is needed to prevent 



zero-frequency modes. However, those theories are lim- 
ited to non-rigid covalent networks or systems interacting 
with attractive potentials, respectively. By contrast, the 
approach adopted here from Ref. [20| can be applied to 
low-coordination crystals, network glasses, as well as re- 
pulsive glasses to explain why all these systems exhibit 
similar structure in the density of vibrational states [13] ■ 

For our model glasses, we have shown that [zi — Zc)/zc 
remains small even when z/zc is arbitrarily large. It is 
for this reason that the marginal-coordination approach 
can be applied to these systems with high coordination. 
As long as (zi — Zc)/zc is small, the behavior of the 
glass is governed by the physics of Point J and isostatic- 
ity. Even though the unjamming transition itself may be 
inaccessible-as it is in the Lennard-Jones system studied 
here-the effect of the longer-ranged part of the potential 
can be treated as a correction. This theory can therefore 
be viewed as a conceptual generalization of the pertur- 
bation theory of liquids to the case of jamming. 
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